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• An algorithm for gauge fixing to the Landau gauge in the fundamental 

^\ ' modular region in lattice QCD is described. The method, a combination 

, of an evolutionary algorithm with a steepest descent method, is able to 

CO ' solve the problem of the nonperturbative gauge fixing. The performance of 

I the combined algorithm is investigated on S**, (3 = 5.7, and IG**, P = 6.0, 

. lattice SU(3) gauge configurations. 
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; 1 Introduction and Motivation 

•i-H . 

' Quantum Chromodynamics (QCD) is the theory that describes the interaction 

^ . between quarks and gluons. From the dynamical point of view, it is usual to 

■ " " ' separate the high energy regime from the low energy regime. While the high 

energy limit of QCD is well described by perturbative methods, perturbation 
theory can not answer a number of important questions. Certainly, it is not 
applicable to the low energy limit of QCD. Presently, we do not have yet an 
analytical method to tackle this dynamical regime. The solution is to solve 
QCD on the computer (1), where continuum euclidean space-time is replaced 
by a discrete set of points, the lattice. Typical lattices are hypercubes where 
points are separated by a in each direction. 

In the lattice formulation of QCD, the gluon fields A° are replaced by the 
links, defined as 

Ufj,{x) = exp{iagoAf,{x + aef,/2)) , (1) 
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where are unit vectors along /i direction. The links are elements of the SU(3) 
group. 

QCD is a gauge theory, therefore the fields related by gauge transformations 

U^,{x) g{x)U^{x) g\x + ae^), geSU{3), (2) 

are physically equivalent. The set of fields related by gauge transformations 
defines a gauge orbit. From the definition of gauge orbit, it follows that to 
study such type of theories it is enough to pick one field from each of the orbits. 
The identification of one field in each gauge orbit is called gauge fixing. 

On the continuum, the problem of the quantisation of gauge theories was 
solved long ago by Feynmanf2), DeWitt (0) and Faddeev and Popov Q). The 
quantisation method requires a choice of a gauge condition, uniquely satisfied 
in each gauge orbit, to define the generating functional for the Green's func- 
tions. For the Landau or the Coulomb gauge and for small field amplitudes, 
the gauge condition is uniquely satisfied in each gauge orbit. However, if large 
field amplitudes are involved, the gauge fixing condition has multiple solutions 
in each gauge orbit (0; 0), the Gribov copies, i.e. the nonperturbative quan- 
tisation of Yang-Mills theories can not be described by the usual methods of 
perturbation theory. This result due to Gribov for the Coulomb and Landau 
gauge was generalized by Singer. In 10), Singer proves that it is impossible to 
find a local continuous and unambigously gauge fixing condition for any SU(N) 
gauge theory defined on the manifold S4. Singer's theorem was extended to the 
four- torus by Killingback Q. 

For the continuum formulation of QCD, in (0) it was argued that the Landau 
gauge Faddeev-Popov formula (5 (9^1) det[— exp[— S'yM(^)], restricted 
to the region where the Faddeev-Popov operator is positive —d ■ D{A) > 
(Gribov region), provides an exact non-perturbative quantisation for QCD. 

The lattice formulation of gauge theories does not require gauge fixing. How- 
ever, gauge fixing is necessary to study the Green's functions of the fundamental 
fields like, for example, the gluon and quark propagators and the quark-gluon 
vertex. The propagators contain information about the mechanisms of confine- 
ment ifiol) and chiral symmetry breaking 0) . The quark-gluon vertex allows 
a first principles determination of the running coupling constant of QCD l|l3t) . 
In addition, by choosing a gauge one can compute renormalisation constants for 
composite operators by sandwiching the operators between quark states 
At least for the usual gauges like the Landau and Coulomb gauges, the lattice 
studies that rely on a gauge fixing condition have a fundamental problem: how 
to properly define a nonperturbative gauge fixing condition, i.e. how to elimi- 
nate the influence of the different Gribov copies on the results. The implications 
of (0) for the lattice formulation of QCD remain to be investigated. 

In this paper we consider the problem of gauge fixing for the Landau gauge. 
On the lattice. Landau gauge fixing can be viewed as a global optimization 
problem |l5). Typically, we have a minimizing function with many local min- 
ima, the Gribov copies, and, to eliminate the ambiguities related to the various 
minima, we aim to find the absolute minimum. In this work, the gauge defined 
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by the absolute minimum of the minimizing function is named minimal Landau 
gauge. 

We report on an algorithm that combines a local optimization method^, a 
Fourier accelerated steepest descent with an implementation of an evolu- 
tionary algorithm ifl^ . suitable for global optimization problems^, to address 
the question of gauge fixing in the minimal Landau gauge. Our investigation 
shows that a proper combination of local and global methods identifies the global 
minimum of the optimizing function and, in this way, solves the problem of the 
nonperturbative gauge fixing. This paper is a full report of the work started in 

ill)- 

The paper is organised as follows. On section 2, the minimal Landau gauge 
is defined. In section 3 the local optimization method, the global optimization 
method and the combined local+global method are described. Section 4 re- 
ports on the performance of combined method for 8^ and 16"* lattices. Finally, 
conclusions and discussion are given in section 5. 



2 The Minimal Landau Gauge 

On the continuum, the Landau gauge is defined by 

a^^A. = 0- (3) 

This condition defines the hyperplane of transverse configurations 

T = {A: d-A ^ 0} . (4) 

It is well known (0) that F includes more than one configuration from each 
gauge orbit. In order to try to solve the problem of the nonperturbative gauge 
fixing, Gribov suggested the use of additional conditions, namely the restriction 
of physical configurational space to the region 

n = {A: d-A = Q, M[A] > 0} C F , (5) 

where A4[A] = —V • D[A] is the Faddeev-Popov operator. However, is 
not free of Gribov copies and does not provide a proper definition of physical 
configurations. 

A suitable definition of the physical configurational space is given by the 
fundamental modular region A C f2, the set of the absolute minima of the 
functional 

Fa[9] = fd^xJ2 Ti:[Af,{x)Af,{x)] . (6) 

The fundamental modular region A is a convex manifold and each gauge 
orbit intersects the interior of A only once (23; 24), i.e. its interior consists 

^By local optimization method we mean an algorithm that seeks only a local solution, i.e. 
a point at which the function is smaller than all other points in its vicinity. 

■^By global optimization we understand the problem of computing the absolute mini- 
mum/maximum of a given function. 
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of non-degenerate absolute minima. On the boundary OA there are degenerate 
absolute minima, i.e. different boundary points are Gribov copies of each other 
I^HSIIB)' The interior of A, the region of absolute minima of identifies a 
region free of Gribov copies. To this choice of gauge we call the minimal Landau 
gauge. 

On the lattice, the situation is similar to the continuum theory I|27tl28l : l29|) . 
The interior of A consists of non-degenerate absolute minima of the lattice 
version of © and Gribov copies can occur at the boundary dA. For a finite 
lattice, the boundary dA, where degenerate minima may occur, has zero measure 
and the presence of these minima can be ignored 

On the lattice, the Landau gauge is defined by maximising the functional 



Fu[g] = ^ Rc { Tr [gix)U^{x)g\x + fi)] } 



(7) 



where 



1 



(8) 



NdunNcV 

is a normalization constant, Ndim is the dimension of space-time, Nc is the 
dimension of the gauge group and V represents the lattice volume. Let (7^ be 
the configuration that maximises F[g] on a given gauge orbit. For configurations 
near on its gauge orbit, we have 



F,7[l + w F,7[l] -f ^ iw°(x)Tr 



X^iU^ix) - U^{x~fL)) - 
A" (C/t(x) - C/t(a:-A))l ,(9) 



where A" are the Gell-Mann matrices. By definition, C/^ is a stationary point of 
F, therefore 



OF 



duj°-{x) 



5]Tr[ A'^(;7^(a;) - [/^(x-A)) 



= 0. 



In terms of the gluon field, this condition reads 

^Tr[ A'^ (A^ix + afi/2) ~ A^ix - afi/2)) 



(10) 



(11) 



or 



0{a) = 0, 



(12) 



i.e. Hl()|l is the lattice equivalent of the continuum Landau gauge condition. The 
lattice Faddeev-Popov operator M{U) is given by the second derivative of (0. 
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Similarly to the continuum theory, on the lattice one defines the region of 
stationary points of Q 

r = {U: d- A{U) = 0} , (13) 

the Gribov's region fl of the maxima of {T)), 

n = {U : d- A{U) = and M{U) > 0} (14) 

and the fundamental modular region A defined as the set of the absolute maxima 
of ((Tj). The lattice minimal Landau gauge chooses from each gauge orbit, the 
configuration belonging to the interior of A. 

The evidence for lattice Gribov copies, i.e. different maxima of Fjj, was 
established long time ago l(30l : Isil but their influence on physical observables 
is not clear. For the lattice Landau gauge, SU(2) simulations suggest that 
the influence of Gribov copies is at the level of the simulation statistical error 
l|33 Ib^ . For SU(3) there is no systematic study but it is believed that the 
Gribov noise is contained within the statistical error of the Monte Carlo. Here, 
we will not discuss the role of Gribov noise on correlator functions but an 
algorithm for finding the absolute maximum of For a discussion on the 

influence of Gribov copies on the gluon propagator see (js^lH^H^. 

3 Optimization Methods 

The algorithm for minimal Landau gauge fixing reported in this paper combines 
a local and a global optimization method. For completeness, in this section we 
outline the local method and describe the global and combined local+global 
algorithms. 

On the gauge fixing process, the quality of the gauge fixing is measured by 
^ = ETr[A(x)At(x)] (15) 

where 

A(x) = ^ [U^{x~ae^) - Ul{x) - h.c. - trace] (16) 
is the lattice version of d^^A^^ — 0. 
3.1 Local Optimization 

By definition, a local optimization method computes a local maximum of ^[/[f;]. 
For Landau gauge fixing, a popular local optimization method is the steepest 
descent (0) method. 

The naive steepest descent method faces the problem of critical slowing down 
when applied to large lattices. Critical slowing down can be reduced by Fourier 
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acceleration. In the Fourier accelerated method, in each iteration one chooses 



g{x) = exp 
where 



p-i i^riB^F VA_. [U^ix) - Ulix)] - trace 



(17) 



A-,((7^(a;)) - C/^(x-ae.) - C/^(x) , (18) 

are the eigenvalues of (—9^), a is the lattice spacing and F represents a fast 
Fourier transform (FFT). For the parameter a we use the value 0.08 (id). For 
numerical purposes, it is enough to expand to first order the exponential in (jl7|l . 
followed by a reunitarization of g{x). 

For large lattices (I17|l is not the best way to solve the problem of critical 
slowing down. In (38; 39) a method was developed that avoids the use of FFT, 
has a dynamical critical exponent close to zero and the advantage to be easily 
parallelized. In this work we use the Fourier accelerated steepest descent method 
(SD). 



3.2 Global Optimization 

Global optimization methods aim to find the absolute maximum or minimum 
of a multidimensional function. Presently, there is not a method that can as- 
sure, with certainty, that the computed maximum in a single run is the absolute 
maximum. Simulated annealing (SA) is, probably, the most popular method for 
global optimization. However, evolutionary algorithms (EA) |l3) are an alter- 
native to simulated annealing. The "advantage" of evolutionary algorithms rel- 
atively to SA is that EA work with multiple candidates for maximum/minimum 
in a single run and, in principle, can avoid or reduce the number of multiple runs 
necessary to identify the global optimum. For us, this provided the motivation 
to try the use of EA for gauge fixing in lattice QCD. 

Evolutionary algorithms (EA) are a generalization of genetic algorithms 
(GA). Genetic algorithms are inspired in natural selection and in the theory 
of evolution of species. The language spoken in evolutionary programming is 
borrowed from genetics. The vector of the parameters to optimize is called 
chromosome or individual. A population consists in a number of individuals. 
The function to optimize is the cost function. 

For the gauge fixing problem, a chromosome is the set of matrices g{x) that 
defined a gauge transformation. The cost function is the functional Fjj [g] . 

An evolutionary code starts generating a set of tentatives of solutions, the 
initial population. In the following the number of individuals in the initial popu- 
lation will be refered by Nipop. In our case, the initial population was generated 
randomly. After sorting the initial population according to their cost function 
value, Npop members were selected, using a rouUette- wheel method (I171') . to be- 
gin the evolutionary phase. The number of individuals in the population was 
always kept fixed to Npop. In this work we used Nipop /Npop = 2.5. 

The population evolution was performed according to the rules 
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1. the best Ngood individuals survive for the next generation; 

2- Nbad = Npop - Ngood are replaced by new chromosomes. The new individ- 
uals are generated by reproducing the Ngood members of the population. 
In this work we set Ngood — Npop /2. 

3. For mating or reproduction two good chromosomes are selected and give 
"birth" to two offsprings. The process completes after generation of 
Nbad new offsprings. In this work parents were selected using the so-called 
roulette- wheel selection (|17f) . a method which favours the best chromo- 
somes in the population. 

4. A new generation is defined only after mutating the population cons- 
tructed after point 3. No mutation was applied to the best member of 
the population. 

In order to reproduce, a population requires a set of rules to make childs 
from the parents, the genetic operators. In this work we considered the following 
operators 

• Random Crossover (RC) 

For evolutionary and genetic algorithms, crossover is a fundamental mat- 
ing operator that mimics the crossover observed in biological systems: 
after selection of a set of contiguous genes in the chromosomes, the two 
childs are built by interchanging the chosen piece of the parents genetic 
material. Our implementation of the crossover is slightly different. On 
the lattice we select randomly V x Pintcross points; the random variable 
Pintcross takcs valucs in [0.40, 0.70]. The offsprings are defined by inter- 
changing the parents matrices g at these points'^ . Note that crossover does 
not implies creation of new genes. 

• Random Blending (RB) 

Blending operators try to overcome the crossover problem of gene cre- 
ation. Our implementation of blending starts by choosing a set of lattice 
points similarly as in RC. For the first child, we select a random value''' 
for (3 G [0,2]. For the selected points, the g matrices are given, after 
reunitarization, by 

Pgi + (1 - /?)ff2, (19) 

where gi and g2 denote parents. In the remaining lattice points we set 
g = gi- The genetic material of the second offspring is generated in the 
same way. The difference being that /? is choosen different at each of the 
selected lattice points and in the remaining points we set g = g2. 

^In literature this type of crossover is also known as uniform crossover. 
*The choice /9 £ [0, 1] is the most simple blending method, but has the disadvantage that 
it does not create new values outside the interval defined by the parents genes. 
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Each mating operator has an associated probabihty. After parent selection, it 
is tested if the chosen parent is able to reproduce by comparing an uniformly 
distributed random number in [0, 1] with the mating probability. 

The mating operators just recombinc the genetic information of the popu- 
lation. To explore more effectively the cost surface, an evolutionary code apply 
mutation operators after the mating phase. These operators change a few genes 
of a chromosome cither by replacing the gene by a neighbouring value or replac- 
ing the gene by a completely different value. In this work, we considered the 
following mutation operators 

• Addition mutations (MA) 

g{x) g{x) + €A, (20) 

• Substitution mutations (MS) 

g{x) A, (21) 

• Expansion mutations (ME) 

g{x) g{x) (1 + eA) (22) 

where e (|e| < 0.025) is a random number and ^ is a random SU(N) matrix. The 
resulting matrix is properly reunitarized. Each mutation operator is applied to 
all population skipping the best A^eHte individuals (in our work Nf,ute = 1). 
Like for the mating operators, mutation have an associated probability too. For 
each operator, and for each individual, we go through the lattice and apply the 
operator in the corresponding g matrices according to the respective probability. 

Each complete iteration of the algorithm (selection, mating and mutation) 
is called generation. 

The probabilities associated to each genetic operator were defined to ma- 
ximize the performance of the pure evolutionary code. The large number of 
parameters makes a detailed study of the probabilities quite hard to perform. 
However, for practical purposes, we used the procedure described below to define 
our algorithm. Set all probabilities to zero except for MS. For MS take 0.01 
for the probability. Change the probability of RC and choose the value that 
optimizes the performance of the algorithm. After setting the probability for 
RC, repeat the procedure for MA, then for RB, then for ME. Finally, check for 
the value of the MS probability. The probabilities associated to the genetic and 
mutation operators used in our study are 
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In what concerns Landau gauge fixing, the performance of the pure evo- 
lutionary algorithm for Landau gauge fixing was quite disapointing; the best 
run ended with 9 ~ 10~^. This can be understood as a consequence of the 
large dimension of the problem and of the nature of Fjj. In conclusion, the 
performance of the pure evolutionary algorithm for the 4D Landau gauge fix- 
ing problem is similar to the performance observed in simplified versions of the 
problem ((ioH ilh. 

3.3 Combined Global + Local Optimization 

For minimal Landau gauge fixing in lattice QCD, the global optimization pro- 
blem can be overcomed by combining the local and the global algorithms de- 
scribed above. From the point of view of the evolutionary algorithm, a possible 
combined algorithm means redefining the cost function as 

f[g;N] = F[/ [5] after local steepest descent steps. (23) 

As described below, with a proper choice of it is possible to identify the global 
maximum of Fjj . 

4 Results for Combined Algorithm and the Min- 
imal Landau Gauge 

The combined algorithm was studied with SU(3) gauge configurations on 8^ 
{(3 — 5.7) and 16^ (/3 = 6.0) lattices. The gauge configurations were gener- 
ated with the MILC code 1^3) using a combination of four over-relaxed and one 
Cabibbo-Marinari updates, with a separation between configurations of 3000 
combined updates. For each of the lattices, the combined algorithm was inves- 
tigated in detail for at least three configurations. 

For Landau gauge fixing, the absolute maximum of Fjj was computed by 
running, for each gauge configuration, 1000 local algorithms for the smaller lat- 
tice and 500 on the larger lattice, starting from different random chosen points. 
A local minimum was defined by demanding that 9 < 10^^" for the smaller lat- 
tice and 9 < lO^^'^ for the larger lattice. The candidate for absolute maximum 
computed with the combined algorithm was compared with the candidate for ab- 
solute maximum from the multiple local algorithm runs. In all the simulations, 
we never observed a larger maximum than the one obtained with the multiple 
runs of the steepest descent method. Preliminary results on the performance of 
the combined algorithm were given in l)2lj) . 

4.1 8^ lattices 

For the S'* lattice, 10 gauge configurations were generated. The study of their 
Gribov copies structure was performed by running 1000 SD on each of the 
configurations. Then, a detailed study of the three configurations with the 
largest number of maxima was performed as described below. 
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Figure 1: Local maxima of one of the 8* SU(3), /3 = 5.7, configurations after 
1000 Steepest Descent starting from random points {6 < 10^^"). 



The number of local maxima computed in the multiple runs with the steepest 
descent method was quite large. Figure ^ resumes the 1000 SD runs for one 
of the configurations used to test the algorithm. The figure shows not only 
a large number of local maxima, the Gribov copies, but also that the most 
probable maxima are associated with the largest values of Fjj. Note that, for 
the configuration considered, the copy with the largest frequency is not the 
absolute maximum. These properties are a general trend observed for some of 
the configurations. 

The combined evolutionary algorithm steepest descent method (CEASD) has 
a large number of parameters and to establish the algorithm we tried to cover, 
as much as possible, the space of parameters. Table ^ is a summary of the 
various runs. All results reported in this paper, for this smaller lattice, consider 
runs with 400 generations and use Nipop ~ 2.5 x Npop , Ngood — Npop /2 and 

Nelite — 1- 

For the combined algorithm, we observed that by increasing N in (|23|l . the 
computed maximum, i.e. the maximum computed after applying a SD to the 
best member of population of the last generation, becomes closer to the absolute 
maximum. Moreover, there is a minimum number of N, Nsteps , such that the 
computed maximum of CEASD is the absolute maximum of Fjj . Figure |21 
reports the number of successful runs of the combined method, for the three 8^ 
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Table 1: Evolutionary populations considered on the 8^ study. The number of 
generations used in each run was 400. 



test configurations, as function of N-ip^p and N . 

Figure 121 shows that, for each population size, there is a minimum value of 
N, Nsteps , such that the CEASD algorithm correctly computes the absolute 
maximum. Figure |21 shows Nsteps as a function of the initial population. The 
solid line is Nsteps for 400 generations and the dashed line is the value of N 
required to identify the absolute maximum in 50 generations. Results seem to 
suggest that for larger populations, Ngteps should become smaller. 

Table |21 reports the first generation that includes the absolute maximum 
in the population^. The results seems to suggest that, for an 8** lattice, 50 
generations may be a safe number of generations for the CEASD algorithm to 
compute the absolute maximum of Fjj . 

Figure reports, for different population sizes, typical evolutions of 9 for 
one of the tested configurations. They show that, in each run, 6 decreases 
rapidly in the first generations, with its value decreasing by roughly 3 to 4 orders 
of magnitude in the first 50 generations, and then remaining approximately 
constant. Moreover, in order to properly identify the absolute maximum of Fu, 
the algorithm seems to require 9 ~ 10~^ — 10^^ after generation 50. 



In conclusion, for an 8"^ lattice it is possible to define a set of parameters 
such that the CEASD algorithm identifies the gauge transformation that max- 
imizes Fjj- For this smaller lattice, our choice being Nsteps — 100, Nipop = 
10 for runs with 200 generations. Note that from figure |21 one reads Nsteps = 
50. However, since evolutionary algorithms are statistical algorithms and the 
combined algorithm requires a relatively low value for 9 to access the absolute 
maximum of Fjj, our choice for Nsteps and the number of generations was con- 
servative. Indeed, results show that similar results can be obtained for runs 
with only 50 generations®. Decreasing the number of generations implies ei- 
ther increasing Nsteps (increasing the computational cost of the cost function), 
increasing N^p^p (increasing the memory requirements) or relying on multiple 
runs of the algorithm. Of course, the user should choose between the different 
possible solutions depending on the computational power he has available. 

^We monitored the presence of the absolute maximum each 50 generations. 

^We tested running the code on 10 configurations for Nsteps = 80, Nipop = 10 and for 200 
generations. Of all the configurations, only one didn't arrive to the absolute maximum. For 
Nsteps = 100, of the 10 configurations tested nine got the absolute maximum in 50 generations 
and only one required 100 generations to compute correctly the maximum. 
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Figure 2: Number of successful runs for the combined method, for the three 8'^ 
test configurations, as function of Nipop and N . 



In order to get an idea on the cpu time required by the CEASD, we bench- 
marked the code on a Pentium IV at 2.40 GHz. For the 8^ lattice, Ngteps = 
100, Nipop — 10 and requiring < 10^^''' for the final steepest descent applied 
to the best member of the population of the last generation, we measured 

number of generations cpu time (s) 
50 2633 
200 10859 

meaning that the CEASD algorithm requires about 54 seconds/generation. For 
the same gauge fixing precision, the steepest descent method requires 56 seconds. 
Therefore, the time required by a run with 200 generations is similar to the time 
required by 200 multiple steepest descent. At this point, a warning should be 
given to the reader. The CEASD code has space for optimization, therefore the 
cpu times reported above should be read as order of magnitudes. The CEASD 
memory requirements for the evolutive phase {Npop — 4) are about 15MB. 
In the next section we report on the CEASD algorithm for a larger lattice. 

4.2 16^ lattices 

For the larger lattice considered in this work, seven /3 — 6.0 gauge configura- 
tions were generated. Similarly to what was done for the 8^ lattice, the Gribov 
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Figure 3: Nsteps versus Nipop for 8"* SU(3), /3 — 5.7, configurations. The solid 
lines gives Nsteps for 400 generations. The dashed line is Nsteps for 50 genera- 
tions. 



copies structure was studied applying 500 steepest descents started from dif- 
ferent randomly chosen points. In order to test the algorithm we performed a 
detailed study for the three configurations with the largest number of Gribov 
copies. 

The first observation being that the number of local maxima is now much 
larger than in the S'' lattice. Figure El shows the Gribov copies found in 500 
multiple SD of one of the configurations used in the detailed study of the CEASD 
algorithm. A similar figure for a 8"* configuration is figure ^ Not only, the 
number of local maxima increases but also the maxima become closer to each 
other when compared to the smaller lattice. To give an idea of the local maxima 
for the 16^ configurations, in table 13 we list the first five highest values of Fu 
computed after the 500 SD method, the 10*^ and the smaller Fu. From the 
numerical point of view, this difference makes the global optimization problem 
a much harder problem to solve. Concerning the frequency of the local maxima, 
the results for the 16^ and 8"^ lattices are similar. The most probable maxima are 
associated with the largest values of Fu but the copy with the largest frequency 
is not always the absolute maximum of Fu ■ 

For the larger lattice, the investigation of the algorithm did not cover the 
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Table 2: Number of generations required by the CEASD algorithm to identify 
the absolute maximum for an 8'* lattice. 
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Table 3: Fu values after 500 SD - 16"^ lattice. 



same set of parameters as in the study of the smaller lattice. Indeed, due to the 
difference on the size of a gauge configuration, a factor of 2^, we only considered 
the smallest population sizes, namely Nipop = 10, 20. The larger populations 
were avoided because of the large memory requirements. In what concerns the 
number of generations on each run, for the larger lattice we only considered runs 
up to 200 generations and checked for the presence of to the best maximum each 
50 generations. As in the previous section, in all runs we used Nipop — 2.5Npop. 

In table^lwe summarize the performance of the algorithm for the three gauge 
configurations considered. For the smallest population, the algorithm seems to 
identify the absolute maximum in 200 generations for N > 180 (Ngteps = 180). 
For runs up to 50 generations, when Nipop = 10 the algorithm sometimes fails 
the computation of the absolute maximum^. For runs up to 50 generations and 
N > 180, the probability of getting the absolute maximum is Pmax = 0.67 for 
Nipop = 10. The probability of getting a maximum which is not the absolute 

''For the remaining 4 configurations, we verified that for N = 180, 190, 200, Nipop = 10 
and for 50 generations the CEASD method computed correctly the absolute maximum of Fij. 
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1 




150 
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1 
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1 


3 


2 


1 


1 


1 
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1 


3 


6 


1 


2 


1 
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1 


1 


2 


1 
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1 


1 


2 
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1 


1 


1 
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1 


1 


1 


1 
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3 


1 


1 
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3 


1 


1 
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200 


1 


3 


1 


1 


1 


1 


160 


50 


1 


1 


2 


2 


1 


1 




100 


1 


1 


2 


2 


1 


1 




150 


1 


1 


2 


2 


1 


2 




200 


1 


1 


2 


1 


1 


2 


170 


50 


1 


1 


2 


1 


1 


1 




100 


1 


1 


2 


1 


1 


1 




150 


1 


1 


2 


1 


1 


1 




200 


1 


1 


2 


1 


1 


1 


180 


50 


1 


1 


1 


2 


1 


2 




100 


1 


1 


1 


1 


1 


2 




150 


1 


1 


1 


1 


1 


1 




200 


1 


1 


1 


1 


1 


1 


190 


50 


1 


> 9 


3 


1 


1 


1 




100 


1 


1 


1 


1 


1 


1 




150 


1 


1 


1 


1 


1 


1 




200 


1 


1 


1 


1 


1 


1 


200 


50 


1 


7 


1 


1 


1 


2 




100 


1 


2 


1 


1 


1 


2 




150 


1 


1 


1 


1 


1 


1 




200 


1 


1 


1 


1 


1 


1 



Table 4: Maxima computed with the CEASD algorithm for three 16^, (3 = 6.0 
SU(3) configurations. 



maximum in K independent runs is then Pother = 0.33 , a number which goes 
rapidly to zero^ with K. Therefore, it seems reasonable to try the use of smaller 
number of generations, provided multiple independent^ runs of the algorithm 
are done. A possible improvement of the multiple run situation could be a 
parallel version of the CEASD algorithm, with the interchange of chromosomes 
between the essentially independent populations every now and then. For the 
largest population, Nipop = 20, the algorithm identifies the absolute maximum 
when 200 generations are considered for N > 170 {Ngt^ps = 170). For runs up 
to 50 generations, again, the algorithm docs not provide the right maximum. 
Now, Pmax = 0.75 for Nipop = 20 and the situation becomes similar to case 
discussed previously. Once more, multiple runs of the CEASD algorithm should 
be able to identify the absolute maximum of Fu when using 50 generations. 

In conclusion, if for runs up 50 generations only a multiple independent run 
can provide the right answer, when the algorithm uses 200 generations, it is 
possible to define Ngteps '■ 

^Pother = 0.33, 0.11, 0.036, 0.012, 0.004 for K = 1, 2, 3, 4, 5. 

^For runs up to 50 generations, if one considers the results for all the seven configurations 
Pmax = 0.86 and Pother = 0.14^. 
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^ ipop ^ steps 

10 180 



20 



170. 



To close this section we report now on the cpu times. On a Pentium IV at 
2.40 GHz, for a 16"* lattice, Nsteps = 200, Nipop = 10 and for 9 < IQ-^^ the 
cpu time measured required by the CEASD algorithm was 



meaning that the CEASD algorithm requires about 2211 seconds/generation. 
For the same gauge fixing precision, the steepest descent method requires 1826 
seconds. Then, the time required by a run with 200 generations is similar 
to the time required by 240 multiple steepest descent. The CEASD memory 
requirements for the evolutive phase {Np^p = 4) are about 236MB. 

5 Discussion and Conclusions 

In this paper we describe a method for Landau gauge fixing that combines an 
evolutionary algorithm with a local optimization method. The "happy mar- 
riage" between the two algorithms is achieved by redefining the cost function of 
the EA, in such a way that it becomes an approximation for the local maximum 
in the neighborhood of the chromosome. In order to get the global maximum, 
the CEASD algorithm seems to require values for of the order of 10~^ for 8^ 
configurations and 10~^ for 16^ configurations. Note that the CEASD requires 
only a good approximation of Fu in order to be able to compute the global 
optimum. 

The combined algorithm was tested for three different configurations in two 
lattices: a smaller 8^ lattice and a larger 16"^ lattice. For both lattices it was 
possible to identify a set of parameters for the CEASD method such that, in a 
single run, the computed maximum, i.e. the maximum obtained after applying 
the steepest descent method to best member of the population of the last gen- 
eration, was always the global maximum defined from multiple steepest descent 
runs. 

For the smaller lattice the CEASD performed extremely well. Indeed, despite 
the relative large number of local maxima, the algorithm seems to be quite stable 
in identifying the global maximum of Fu - see table|21 For the larger lattice, the 
number of local maxima is much larger when compared to the 8"* lattice. Not 
only the number of maxima is larger but they are closer to each other. From the 
point of view of the global optimization, this means that the numerical problem 
in hands is much harder to solve. Nevertheless, again it was possible to define 
a set of parameters such that the algorithm identified the global maximum 
in all tested configurations - see table 0| Our choice of parameters for the 
CEASD algorithm (200 generations, Nsteps — 100 for 8^ lattice and Nsteps = 
200 for the larger lattice and for Nipop = 10) is a conservative choice. As 



number of generations 



cpu time (s) 



50 
200 



112090 
436071 
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explained before, it is possible to use smaller values of N or smaller number of 
generations. Decreasing N and/or the number of generations implies decreasing 
the run time of the CEASD algorithm. However, reducing N and/or the number 
of generations should be done with care. Indeed, the comparative study of the 
two lattices shows that the complexity of the maximisation problem increases 
with the lattice size, that the method works better for larger populations, larger 
values of N and for sufficiently large number of generations. Nevertheless, the 
results of the previous section also show that, for relatively large values of N, the 
probability of computing the absolute maximum of Fu is large. This suggests 
that a possible solution to the global optimization problem is to perform multiple 
independent CEASD runs using lower values of N and/or smaller number of 
generations. For sufficient number of independent runs, in principle, the method 
should be able to get the global maximum. A similar situation is found when 
one relies on simulated annealing for global optimization problems. 

The cpu times required by CEASD algorithm for the two lattice sizes seems 
to suggest that the scaling law of the combined method is close to the Fourier 
accelerated SD method, i.e. \nV with 6 taking values close to 1. A measure 
of S requires necessarily an analysis with more lattice sizes^". This is a numerical 
intensive problem. We are currently engaged in measuring S and will report the 
result elsewhere. Naively, one expects that gauge fixing to the minimal Landau 
gauge with CEASD is as demanding as performing a gauge fixing with the SD 
method. 

In principle, it is possible to combine the EA with any local optimisation 
method. Faster local methods will produce faster combined algorithms. The 
time required by a combined algorithm is strongly dependent on the perfor- 
mance of the local method. The gauge fixing is a computational intense prob- 
lem. Therefore, it is important to investigate new and more performant local 
methods. 

The CEASD algorithm described here for Landau gauge fixing seems to 
solve the problem of the minimal Landau gauge fixing. Moreover, the method is 
suitable to be used with other gauge conditions that also suffer from the Gribov 
ambiguity and are currently used in lattice gauge theory. The effects of Gribov 
copies in QCD correlation functions remains to be investigated l|36|) . 

P. J. S. acknowledges financial support from the Portuguese FCT. This work 
was in part supported by fellowship Praxis/P/FIS/ 14195/98 under project Op- 
timization in Physics and in part from grant SFRH/BD/10740/2002. 

^"Assuming a scaling law like InV and using the cpu times reported in this work, we 
get (5gj-) = 1.15 for the Fourier accelerated method and ^Q-g^gj-j = 1.24 for the CEASD 
algorithm. For a 32'* configuration, these numbers mean that a SD gauge fixing requires about 
T = 6 X 10'* s, the CEASD requires 72 X r seconds to run in 50 generations and 282 X r seconds 
for a 200 generation run. 
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Figure 4: log(^) for one of the 8^ test configurations and different population 
sizes. The value of N is represented by the thickness of the line (larger thickness 
meaning larger N). 
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Figure 5: Local maxima of a 16* SU(3), /3 = 6.0 configuration achieved after 
500 Steepest Descent starting from random starting points {9 < 10~^^). For the 
three configurations used to set the algorithm, the number of different Gribov 
copies found after 500 SD was 177/500, 238/500, 326/500 for configurations 
number 66000, 72000 and 9000, respectively. 
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